%% OSLWORKSHOPS - DEMONSTRATING HOW TO REMOVE SCANNER ARTEFACTS

% This is an example of how to remove Elekta Scanner Artefacts from
% your data for successful MaxFiltering.

% You will use
% osl_call_maxfilter.m - a wrapper function for calling MaxFilter.
% osl_convert_script.m - a function to convert .fif files into SPM objects.
% oslview.m            - OSL's continuous data viewer.

% You will need OSL version beta.2 or later (earlier versions may work but no
% guarantees).

% and the practical data files

% rest_data_raw.fif
% spm_rest_data_raw_nosss.mat/.dat
% spm_rest_data_raw_sss_nocorr.mat/.dat

% Henry Luckhoo
% henry.luckhoo@trinity.ox.ac.uk
% 13.11.12

%% Setup Paths

osldir='/Users/woolrich/homedir/matlab/osl1.2.beta.3';
rmpath(genpath('/net/horus/data/OHBA/spm8')); % If you have other versions of SPM on your path you should remove them too!
rmpath(genpath('/opt/analysis/spm8')); % If you have other versions of SPM on your path you should remove them too!
addpath(osldir);
global OSLDIR;
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set-up files and directories

workingdir = '/home/hluckhoo/Desktop/ohba_workshop_artefact_rejection_data/osl_example_maxfilter_data'; % set this to the directory with the 5 data files listed above.
cd(workingdir);
fif_files = {[workingdir '/rest_data_raw']};
spm_files = {[workingdir '/spm_rest_data_raw']};

%% Let's look at the raw data using OSLVIEW

D=spm_eeg_load([spm_files{1} '_nosss']);
oslview(D)

% 1.) In the GUI, change the time axis to look at 2 or 3 minutes of data by
% clicking on the "Increase time window button" in the top left (note that
% the x-axis is labelled in seconds). 

% 2.) Have a look at the data. Click on the "P" button to switch to
% Magnetometers "M".

% 3) Can you see the single corrupted magnetometer with a spike and then
% decay. Close the GUI. When prompted to save the figure, select "No".

%% Let's look at some data where we have applied Maxfilter using OSLVIEW

D=spm_eeg_load([spm_files{1} '_sss_nocorr']);
oslview(D)

% 1.) In the GUI, change the time axis to look at 2 or 3 minutes of data by
% clicking on the "Increase time window button" in the top left (note that
% the x-axis is labelled in seconds). 

% 2.) Have a look at the data. Click on the "P" button to switch to
% Magnetometers "M".

% 3.) Can you see how the artefact has been spread to all the channels. 
% Close the GUI. When prompted to save the figure, select "No".

%% How to prevent this - an example pipeline for avoiding these
%% artefacts! Use the following as a model for your own analysis.

% You may need to clean up any previous Maxfilter results (if it has been run
% before). If so run something like this:
% runcmd(['rm ' fif_files{1} '_sss.fif ' fif_files{1} '_nosss.fif']) % This is just to delete any partially processed fif_files.

%% 1 - Create Un-maxfiltered data & Convert to SPM

for subnum=1:length(fif_files) % Loop over subjects
    Smf=[];
    Smf.fif=[fif_files{subnum}];
    Smf.logfile=1;
    Smf.downsample_factor=4;
    Smf.nosss=1;
    osl_call_maxfilter(Smf);
end


for subnum=1:length(fif_files), % Loop over subjects
    S=[];
    S.fif_file=[fif_files{subnum} '_nosss.fif'];
    S.spm_file=spm_files{subnum};
    osl_convert_script(S);
    close;
end

%% 2 - Mark any channels as bad in the SPM object using OSLVIEW

D=spm_eeg_load(spm_files{1});
oslview(D);

% 1.) In the GUI, change the time axis to look at 2 or 3 minutes of data by
% clicking on the "Increase time window button" in the top left (note that
% the x-axis is labelled in seconds). 

% 2.) There are two noisy planar gradiometers but MaxFilter can deal with
% them. Click on the "P" button to switch to "M" Magnetometers. Can you see
% a channel with a massive jump and then decay. MaxFilter cannot handle
% that. So set that channel to bad by right clicking on the noisy channel 
% and selecting "Set Channel as Bad". 

% 3.) Click on the save buttton and then close the GUI. When prompted to
% save the figure, select "No".

%% 3 - Let's check that we have saved the bad channles successfully

clear D;
D=spm_eeg_load(spm_files{1});
D.badchannels

%% 3 - Re-run Maxfilter BUT use this badchannel information

for subnum=1:length(fif_files)
    Smf=[];
    Smf.fif=fif_files{subnum};
    Smf.logfile=1;
    Smf.downsample_factor=4;
    %Smf.nosss=0;                    % Smf.nosss=0; is the default value.
    Smf.spmfile = spm_files{subnum}; % We pass in the name of the SPM file!
    osl_call_maxfilter(Smf);
end

S=[];
for subnum=1:length(fif_files), % iterates over subjects
    S.fif_file=[fif_files{subnum} '_sss.fif'];
    S.spm_file=[spm_files{subnum} '_sss'];
    [D spm_files_sss{subnum}] = osl_convert_script(S);
    close
end

%% 4 - Check the resulting output

% Raw Data
Draw=spm_eeg_load([spm_files{1} '_nosss']);
oslview(Draw);

% Maxfiltered Data without flagging bad channels
Dnocorr=spm_eeg_load([spm_files{1} '_sss_nocorr']);
oslview(Dnocorr);

% Maxfiltered Data with bad channels flagged
Dcorr=spm_eeg_load([spm_files_sss{1}]);
oslview(Dcorr);

% See what a difference getting rid of that bad channel made. 